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Abstract 

Numerical models are starting to be used for determining the future behaviour of 
seismic faults and fault networks. Their final goal would be to forecast future large 
earthquakes. In order to use them for this task, it is necessary to synchronize each 
model with the current status of the actual fault or fault network it simulates (just 
as, for example, meteorologists synchronize their models with the atmosphere by 
incorporating current atmospheric data in them). However, lithospheric dynamics 
is largely unobservable: important parameters cannot (or can rarely) be measured 
in Nature. Earthquakes, though, provide indirect but measurable clues of the stress 
and strain status in the lithosphere, which should be helpful for the synchronization 
of the models. 

The rupture area is one of the measurable parameters of earthquakes. Here we 
explore how it can be used to at least synchronize fault models between themselves 
and forecast synthetic earthquakes. Our purpose here is to forecast synthetic earth- 
quakes in a simple but stochastic (random) fault model. By imposing the rupture 
area of the synthetic earthquakes of this model on other models, the latter become 
partially synchronized with the first one. We use these partially synchronized mod- 
els to successfully forecast most of the largest earthquakes generated by the first 
model. This forecasting strategy outperforms others that only take into account 
the earthquake series. Our results suggest that probably a good way to synchronize 
more detailed models with real faults is to force them to reproduce the sequence of 
previous earthquake ruptures on the faults. This hypothesis could be tested in the 
future with more detailed models and actual seismic data. 

Key words: Earthquake prediction, fault model, cellular automata, 
synthetic-earthquake catalogues, seismic modelling, characteristic earthquakes. 
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1 Introduction: Data assimilation in dynamical fault models 



Numerical mo dels are now frequently used to simulate the seismic behaviour 



of fa ults (e.g. Kato and Send. 2003c Fitzenz and Millei . 2004 : Kuroki et al 
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2004 ) and fault networks ( e .g.lWardi 1200 0: Hashimoto, 2001; Robi nson and Benitesl . 



200l[ iR.undle et all 1200 it ISoloviev and Ismail- Zadeb . ■2003t .Robinsonl . l20ol 



Rundle et al.l . l2004^ . In these models, fault planes separate lithospheric blocks 



that are strained at specific rates, and sudden slips (earthquakes) are gener- 
ated by the faults according to certain friction and/or rupture laws. Although 
no completely realistic dynamical model presently exists, these simulations 
are now sufficiently credible to begin to play a substanti al role in scientific 
studies of earthquake probability and hazard ( Ward , 2000l ) . The final goals of 
the numerical modelling of seismicity are not different from, for example, the 
goals of numerical models of the atmosphere. A good model should be able 
to: 



(1) reproduce the general characteristics of the system, 

(2) mimic the state of the system at the present moment, and 

(3) forecast the future evolution of the system. 

Most numerical models of seismicity have been designed to achieve the first 
goal, by reproducing gener al characteristics of earthquakes such as their size- 
frequency distribution (e.g . Bak and Tand. 19891: 



19981 : IPreston et aL 



200 



Olami et al.lll992l : iDahmen et al 



Vazquez-Prada et al 



2OO2I ). or the generation of 



aftershock and foreshocks (e.g. iHainzl et al. . ll999l ). When a model is designed 
this way, it is left to evolve freely according to its rules, and all that is checked 
is whether the overall results of the model are similar to the observations made 
in Nature or not. 



The second goal requires data assimilation, that is, the process of absorbing 
and incorporating observed information into the model. By this process, the 
model is tuned and synchronized, at least partially, with the real system it tries 
to simulate. In a meteorological model, data of atmospheric pressure, temper- 
ature, humidity, cloud cover, precipitation, etc. measured in a given moment 
at different locations and heights can be included. With this procedure, the 
model becomes a reasonably good representation of the atmosphere at that 
moment. Then it can be used to calculate the probable future atmospheric 
evolution (i.e. the third goal cited above). 

Seismic data assimilation poses greater problems than its meteorological equiv- 
alent. This explains (at least partially) the relative delay in developing reli- 
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able for ecasts of large e arthquakes. The inner workings of both the atmo- 
sphere (i Houghtonl 120021 ) and the lithosphere (|GoltzL Il997l : iTurcotteL Il997 : 
Keilis-Borok. .2002) are complex and chaotic, so they are inherently difficult to 
forecast. However, while meteorologists can probe the atmosphere every day at 
different places and heights (and assimilate the obtained data in their models 
in near real-time), lithospheric variables of paramount importance, such as the 
stress and strain, can be measured only in certain places, and n ot at any time: 



earthquakes have unobservable dynamics ( Rundle et al. . 2003f ). For example 



the best current compendium of stress naagnitudes and directions in th e litho- 
sphere is the World Stress Map ( Zoback , 19921 : Reinecker et al. , 2004 ) , whose 
entries are point static time-averaged estimates of maximum and minimum 
principal stresses in space. And the dire ct measurements of stress on active 
fault zones at depth are still sc arce fe .g. Ilkeda et al.l I2OOI"; 'Tsukahara et al 



2001 



Yamamoto and Yabd 12001: Hickman and Zoback . .2004: Boness and Zoback , 



2004^ . The dynamical models would need better spatial and temporal infor- 
mation of stress, both more abundant and m ore systematically collected than 
that currently available (|Rundle et al.Ll2004^ . It is thus necessary to seek ways 



to tune and synchronize the models with more abundant observable data. 



A first step of data assimilation in models of earthquake faults is to intro- 
duce information regarding the topology (that is, the shape and location) of 
the active faults and their long-term behaviour. For example, the long-term 
fault slip rate, and the average recurrence interval of the largest earthquakes 
in the fault can be estir nated from paleos e ismol ogical studies and should be 
included in the models ([Grant and Gouldj . 200^ ). Examp l es of this approach 
are the works of iR.undle et all (|200lL bool and iRobinsonI (|2004l) . The surface 



deformation measured via Global Positioning System (GPS) networks and by 
Synthetic Aperture Radar Interfero nietry (InSAR ) can also constitute input 
data for the dynamical fault models ( Rundle et al.L 120041 ). Earthquakes them- 
selves are indeed the most obvious observable events of lithospheric dynamics, 
and could provide the most detailed data available to assimilate in the models, 
but how? The earthquake rupture area could be an important clue. 



The r uptu r e area and slip distribution in real e arthquakes can be very com- 
plex ( Sieh . 19961 : Kanamori and Brodskv . 2004), but can be estimated in a 
variety of ways. The actual s lip distribution can be obtaine d by inverting the 
observed s eismic waveform s ("Kana mori and Brodskv , 2004 ) or tsunami wave- 
forms (e.g. Tanioka et al. . 2004iBaba and Cummins . 2005[). and/or by geode- 
tic modeUing of surface displacement ( Yabuki and MatsuuraL Il992[ ). Some 
earthquakes produce surface rupt ures, which are useful for estimating the 
rupture area ( Stirling et al. . 2002). Although most surface ruptures occur in 
large shocks, with magnitudes larger than about 6, they have been reported 
for earthquakes with magnitudes dow n to 2.5 (see the co mpilation of historic 



earthquakes with surface rupture by lYeats et al.L 119971 pp. 473-485). Also 



the rupture area can be estimated from the seismic moment (calculated from 



3 



the a mplitude spectra of seismic waves- IScholzl. 120021: iKanamori and Brodskvl. 



2004 ) , or from the moment m a gnitu de ( Wells and Coppersmithl . ll994 ; 



Stirling et al. 



2002; Dowrick and Rhoadesl . 2004). Frequently the locatio n of early after- 



shock s is used to determine the rupture area of the mainshock flWells and Coppersmith , 
1994 ). although the aftershock zone tends to grow wi th time (iKisslingerl 19961 ) 



and is not necessarily a good indicator of that area ( Yagi et all 1999| ) 



Complex models with realistic fault topology are able to reproduce the rup- 
ture area and coseismic slip of historical earthquakes. It is thus possible to 
force the model to reproduce the rupture of a historical earthquake, and let it 
evolve fr om tha t morn ent onwards to see what could happen in the future. For 
example, Ward ( 2000l ) developed a model including the network of main faults 
in the San Francisco Bay Area (California). He forced the model to reproduce 
the San Andreas Fault surface coseismic slip of the 1906 San Francisco earth- 
quake, and let it evolve freely from that earthquake onwards, in an attempt to 
simulate the probable sequence of earthquake ruptures during the next 3000 
years. 



But considering only the data of the largest earthquake in the series is probably 
not sufficient to properly synchronize the model. Complex and chaotic systems 
are very sensitive to the initial conditions. The information regarding only one 
event probably does not sufficiently constrain the initial conditions, and the 
calculated evolution will probably be a particular case of a large range of 
possible outcomes. Will this panorama improve by forcing the model to repro- 
duce all the observed earthquake ruptures, including the small ones? Probably 
yes. To check whether this idea works, at least to forecast synthetic seismic- 
ity, is the purpose of this paper. The number of recorded large earthquakes 
is relatively scarce, especially in individual faults (where the recorded series 
very rarely includes ten large events). This hampers the ability to characterize 
statistically the effectiveness of any forecasting method. Synthetic earthquake 
catalogues, on the other hand, can be as long as desired. This enables to as- 
certain, with robust statistics, whether a forecasting strategy could be useful, 
before endeavouring to apply it to real seismicity. 



In the following sections, our goal will be to forecast the largest earthquakes 
generated by the minimalist model, a simple numerical fault model. We will 
show that when all the earthquake ruptures generated by this model are im- 
posed on other, similar models, these become partially synchronized with the 
former. We use them to declare alarms that efficiently mark the occurrence of 
the largest shocks in the first model. The results are much better than those 
obtained with other strategies that consider only the earthquake series. The 
model, albeit simple, is stochastic (it involves randomness), so its efficient 
forecasting is not trivial. We will describe how this stochasticity can be dealt 
with, by using an approach similar to the so-called ensemble forecasting used 
in Meteorology ( Palmer et al. . 20051 ). The method could be used in other more 
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detailed and realistic models (stochastic or not) to test our general conclusion: 
that they might be partially synchronized with actual faults by being forced 
to reproduce the series of observed earthquake ruptures. 

In the next section we describe the model and its properties. Then, we outline 
the general scheme of prediction and the forecasting strategies used as refer- 
ence to assess the merits of any other predictive method in the model. Finally, 
the method based on partial synchronization is explained and its possible 
utility discussed. 



2 The minimalist model 



The minimalist model is the numerical model whose l argest earthquakes we 



will t ry to forecast. It was introduced in a previous work (jVazquez-Prada et al 



2002| ). and has mainly two, apparently contradictory, advantages for the pur- 



pose of this paper: it is simple but, at the same time, it is difficult to forecast. 
Because it is simple, several of its properties can be derived analytically, and 
it can be characterized in detail with numerical simulations which do not re- 
quire an impractical amount of computer time. Because it is stochastic, it is 
difficult to forecast, so the results we will obtain here are not trivial. In the 
following paragraphs we will explain how the model works, and what are its 
main properties, comparing them with those of actual faults. 



2.1 How the model works 



The model is a simple (hence its name) cellular automaton. Cellular automata 
are frequently used to model seismic faults. In these models, the fault plane 
is divided into a grid of cells (each cell representing a fraction of the fault 
area), and the time evolves in discrete time steps. Each cell's state is updated 
at each time step according to rules that usually depend on the state of the 
cell or that of its neighbors in the previou s time step. The se rules can be 
designed according^ to certain fr i ction laws (iBen-Zion j_ 2001 ) , stress transfer 
flOlami et al.Lll9 92: Ha inzl et all llQQfll : IPreston et al.l . I2OOOI) . and the effect of 
fluids f)Miller et al.L 119991 ). In the minimalist model, as well as in other very 

simple cellular .nU^^. (e.g. i New,... ..H T„rnntt,J H I 

iGonzalez et al.L 

200i), these details are ignored: the model is driven stochastically, there are 



only two possible states for each cell, and the earthquakes are generated ac- 
cording to simplified breaking rules. 



Let us now explain the simplified view of earthquake generation that the model 
tries to sketch. In actual faults, the regional stress strains the rock blocks of the 
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fault, making portions (patches) of the fault plane to become metastable. That 
is, they are static, but store enough elastic energy to propagate an earthquake 
rupture once triggered. Different processes (for example, fault creep -aseismic 
slip- and plastic deformation) dissipate stress along the fault plane, so stress 
is not directly converted into elastic strain. Earthquakes rupture some of the 
metastable patches of the fault, that then become stable, thus relieving strain. 
The hypocentre of an earthquake is usual ly located in a particularly s t rong 
patch of the fault plane, called "a sperity" ( Kanamori and Stewart . 19781 : Aki, 
19841: |d3, bOO.i Ihei et ahLliil . Asperities appear to be pers i stent features 



where earthquake ruptures start once and again (jAkil . 119841 lOkada et al. 



2003| ). Once the rupture starts, it propagates along the fault plane until it 



arrives at a patch of the fault that is not sufficiently strained. Then the rup- 
ture cannot propagate further, and is arrested. The relatively stable patch that 
i s not sufficiently strained and that arrest s the rupture is called the "barrier" 
f)Das and Akil . Il977l : IXkil . Il984l : 1513 . 120031 ). 



The model, depicted in Fig. 1, sketches these features as follows. It divides 
the plane of a fault into an array of N equal ce lls, eac h denoted by an index i. 
In previous papers (IVazauez-Prada et al. . 2002i..2003i : Lopez- Ruiz et al. . 2004; 
Gomez and PachecoT 200^ Gonzalez et al.l 12004 ). this array was drawn ver- 



tically, in order to simplify its mathematical description. Here the model will 
be drawn horizontally, in order to sketch the fault plane in a way more sim- 
ilar to that of actual faults (which are usually longer along the strike than 
along the dip). Some oth er cellular automato n models discretize the fault 
plane in a similar way fe.g. iRundle et aL . 2004). The parameter N is the only 
one that can be changed in the model. The cells can only be in one of two 
states: "empty" (stable) or "occupied" (metastable). The state of the model 
at each time step can be described simply by stating which cells are occupied 
and wh ich are not. The increase of regional stress, as in other simple models 
feak and Tane . 19891 : [Newman and Turcottd . 20021 : Castellaro and Mulargia . 
EoOll : iGonzalez et al.l . Eiol T is represented by the random addition of "stress 
particles". This randomness is a way of dealing with the complex stress in- 
crease in actual faults. At each time step, one cell is selected randomly, and 
a new particle arrives on it. That is, each cell has the same probability, 1/A^ 
of receiving the new stress particle. If the chosen cell is empty, the particle 
"occupies" it. This means that the regional stress has produced enough strain 
on that cell to make it metastable. If the cell is already occupied, that stress 
particle is lost; this is analogous to stress dissipation on the fault plane. The 
total number of occupied cells represents the total elastic strain on the fault. 



In the model, we assume that there is only one, persistent, asperity: the first 
cell, i = I, placed at one end of the array. This option is chosen because it 
simplifies the analytical description of the model. When a stress particle fills 
cell i = 1, a. rupture starts there, and propagates through all the consecutive 
metastable cells until it is arrested by a stable cell. That is, if all the successive 
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/= 1, 2 N 



Fig. 1. The minimalist model as a sketch of a seismic fault. The fault plane is divided 
into an array of N equal cells, denoted with an index z, 1 < i < A^. The increase 
of regional stress is represented by the random addition of "stress particles" to the 
cells. Earthquake ruptures start at an asperity, the cell i = 1, when a stress particle 
arrives to it. The rupture propagates through all the consecutive metastable cells 
(occupied by particles). The rupture area is k, the number of cells broken. The 
figure depicts an earthquake with k = 3. 



cells i = 1 to i = k are occupied, and cell + 1 is empty, then the effect of 
the earthquake is to empty all the cells from i = 1 to i = k. The other 
cells, i > k remain unaltered. The cell /c + 1 is a barrier: it is empty (stable), 
so the rupture cannot propagate through it. The size (rupture area) of the 
earthquake is k, the number of cells broken in the synthetic earthquake. Thus, 
the earthquake size in the model is discrete, 1 < A; < A^. Earthquakes, in 
practice, are instantaneous in the model (they do not last for any time step). 
This represents the fact that earthquake ruptures are, indeed, much faster 
than the slow stress loading represented by the addition of particles. 

The random addition of particles is what makes the model stochastic. It also 
determines the rate at which earthquakes occur in the model. At each time 
step, independently of the previous earthquake history, there is a probability 
1/N for the incoming stress particle to arrive at cell i = 1 and start an 
earthquake. Thus an earthquake, on average, occurs every steps. The time 
between any two consecutive earthquakes is purely random (Poissonian, with 
rate 1/A^). 



The cellular-automaton approach of this model is similar to that of the "for- 
est fire" models, in which clusters of interconnected occupied cells ("trees") 
"burn" and are reset to empty when they are randomly struck by "lightning" 
( Drossel and Schwabll . 11992: lHenlevl . E993[ ). The utility of this k ind of models 
for earthquake physics has been noted by Rundle et al. ( 2003[ l. In the min- 
imalist model there is no random "lightning": the clusters of interconnected 
metastable sites are only emptied if they are connected to the cell i = 1 and 
if this fails. 
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2.2 Main properties of the model 



The minimahst model, because of its extreme simplicity, lacks the detailed 
description of the seismic process that a fully dynamical model can display. 
For example, it does not include the effects of fault friction, elastic stress 
transfer, or the role of fluids that more complex models can take into account. 
However, it spontaneously displays several properties that are comparable to 
those of actual faults, outlined as follows: 



(1) Eart hquake size-frequer i cv distribution. It is of the characteristic-earthquake 



(2) 



(3) 



type ( Wesnouskv et al. . l983[ Schwartz and Coppersmith . 1984 : Youngs and Coppersmith . 



Ll985l; 'Wesnouskv', 1994'), observed in seismic faults with simple traces 
(^Stirling: et al.. . 1996) and in other nu r nerical models if the fault plane is 
homogeneous (e.g. iR.undle and Kleinl. Il993l: iMainl. 11996^ 

iDahmen et al 



1998 [lsteacv and McCloskevl.ll999l : lMoreno et al.Lll999l : lHainzl and ZoUer 



2nni[lHeimDer 120031: IZoUer et al.L l2005| ). In this distribution there is a 



relative excess of events (called characteristic earthquakes) which break 
the whole fault or most of it. In the model, they are the earthquakes with 
size iV, a nd will be the events to f orecast. The Gutenberg - Richter dis- 
tribution ( Ishimoto and "Tkkl . I1939I : iGutenberg and Richteil . Il944 Il954l ) 
observed in regional seismicity (which includes contributions from many 
faults) can be reproduced adding up the seismicity of an ensemble of 
minimalist models whos e sizes (A^) are distributed as in actual faults 
(|L6Dez-Ruiz et al.l . l20oJ ). 

Duration of the earthquake cycle. The earthquake cycle of a fault is the 
time interval between t wo consecutive characteristic earthquakes (e.g. 
IScholz and Gupta . 2000| ). The statistical distributio n of these intervals iii 
the model is similar to the observed in seismic faults ( Gomez and Pachecol . 



(4) 



2004^ . The distribution of time intervals between consecutive earthquakes 
of any size in the model is Poissonian distributed. However, if only the 
characteristic earthquakes are considered, the distribution is not Poisso- 
nian. This is because the maximum possible size of an event depends on 
the size of the previous event and the time elapsed since it occurred. This 
is commented in the next paragraph. 

Stress shadow. When a fault generates a large earthquake, the elastic 
strain is reduced, and a minimum time has to elapse until the fault, by 
slow tectonic deformation, accumulates enough strai n to generate another 
large earthquake. This effect is called stress shadow ( HarrisL 2000( ). In the 
minimalist model there is a stress shadow: if an event of size k takes place, 
at least k time steps have to elapse until another event of that size can 
occur. 

Pattern of strain loading. In actual faults, the stra in increases rap idly iust 
after a large earthquake, and then more slowly (|Michael l200,5l ). In the 
model, the total elastic strain is represented by the occupation (the total 
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400 600 800 
Time steps 



1000 1200 



Fig. 2. The number of occupied cells in a minimalist model with = 20, for ten 
seismic cycles. This number is analogous to the total elastic strain accumulated in 
the fault. Sudden drops correspond to earthquakes. Each seismic cycle ends with 
an earthquake of size A^. 

number of occupied cells, Fig. 2), which has a similar pattern. Just after a 
large earthquake, there are fewer occupied cells, so it is more probable for 
the incoming particles to land on empty cells, and the occupation grows 
faster than later on (Fig. 2). 
(5) Seismic quiescence. The model displays seismic quiescence (absence of 
earthquakes) before the characteristic events. Once N — 1 cells (from 
i = 2 to i = N) become occupied, the occupation reaches a plateau (Fig. 
2) and, on average, time steps have to elapse until the next earthquake 
(which is the characteristic one) occurs. Seism ic quiescence has be en 



[e.g. 



Wvss and Habermann 



observed preceding many large earthquakes ^ ^ _^ , 

fl988; .Scholz. . .2002) , although in other cases the opposite effect (increased 
activity) has been ob served (e.g. Bowman et al.l . 19981 : Reasenbere . 19991 : 



Tiampo et al. . 20021 ). The minimalist model does not show this last be- 



haviour. 



3 General scheme of forecasting 



In this section we will explain the general framework for the forecasting of the 
largest earthquakes in the model. As a first remark, we have to consider that 
the model is stochastic, so it is not predictable with absolute precision. Only 
simple deterministic systems are fully predictable. The evolution of complex 
systems, such as the atmosphere or the lithosphere (even if it were determinis- 
tic) is very sensitive to the initial conditions. As these complex systems cannot 
be fully characterized, they turn out not to be fully predictable either. 



Earthquake predict ion ( Keilis-BorokL 20021: Keilis-Borok and Solovie v. 2UU^ 
Rundle et al.l . 2003| ) , as well as some atmospheric predictions (|Masonl . 120031 ) 
is frequently regarded as a binary forecast: one has to decide whether a large 
earthquake is going to occur or not, in a certain time-space window, instead 
of calculating the exact probability of this event. In this binary-forecasting 
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approach, an "alarm" is declared when a large earthquake is expected. If it 
takes place when the alarm is on, the outcome is a successful forecast. If it 
takes place when the alarm is off, there has been a prediction failure. If the 
alarm was declared during a certain period, but the expected earthquake did 
not happen, that constitutes a false alarm. 

Note that for using this approach it is necessary to define precisely what the 
target earthquakes are that we wish to forecast. Usually they are defined as 
those with a magnit ude larger than a given threshold, both when dealing w ith 
actual earthquakes dKeilis-Borok and Soloviey , 20031: Rundle et al.L 200^ or 



Pepke and Carlsonl . ll994l : lHainzl et al.l . l2000h . In the 



with synthetic ones (e.g. 
minimalist model, it is natural to choose as target events the characteristic 
earthquakes (size k = N), as they mark a distinct peak in the size- frequency 
diagram, being much more frequent than other large earthquakes. 

A way to quantify the forecasting ability of a certain strate gy is to comput e 



the fraction of errors, fe, and the fraction of alarm time, fa (iMolchanl . 119971 ). 
Given a certain time series of the model, fe is the ratio of the total number 
of prediction failures to the total number of target events. And fa is the ratio 
of the total time during which the alarm was on to the total duration of the 
time series. The fraction of false alarms, ff, is included in fa, and is the ratio 
of the total duration of false alarms to the total duration of the time series. 

Of course, a good forecasting strategy should render small fa, fe and //. 
However, as a general rule, a strategy that renders low fe tends to produce 
large fa and //. Dealing with real seismicity, both a failure and an alarm are 
costly. Eventually, decision-makers would need to consider what is less costly: 
to predict most of the dangerous earthquakes, but declaring many alarms, or 



to de clare fewer alarms but failing the forecast of more large shocks (jMolchan 



19971). Depending on the trade-off between costs and benefits, one should try 



to minimize a loss function, L, that can depend on fa, fe and/or ff. 

In the next section, we will describe the forecasting strategies that will be used 
to compare the merits of the new strategy proposed in this paper, based on 
synchronizing models between themselves. In the first of the subsections we 
will indicate the loss function we will try to minimize in the forecasting of the 
model. 



4 Forecasting strategies for comparison 



We describe here three forecasting strategies, based on the earthquake series, 
that we will use to asses the merits of the new strategy described later in 
this paper. The first two strategies (the random guessing strategy and the so- 
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called reference strategy) can be used in any system. The third is specific to 
the minimalist model, and serves to ideally determine its maximum theoretical 
predictability. 



4-.1 Random guessing strategy 



In this strategy, the alarm is randomly turned on and off, during a certain 
fraction of alarm time, /„. It is simple to apply this strategy to any cellular 
automaton model. Here, in each time step, the alarm is on with a probability 
p. As a result, the alarm will be on during a fraction p of time steps (/„ = p). 
When the target earthquake finally occurs in a certain time step, there will be 
a probability p for the alarm to be on. Thus, on average, a fraction p of target 
earthquakes will be predicted, and a fraction fe = 1 — p will be prediction 
failures, so /a + /e = 1. This strategy has two trivial cases: if the alarm is 
always on (/„ = 1), all the target earthquakes are "forecasted" (/e = 0). 
Conversely, if the alarm is always off {fa = 0), we fail to predict any of them. 

To be statistically significant, any forecasting strategy must render better 
results than a random guess. A natural way to measure this improvement is 
to consider the loss function L = fa + fe- Then, L = 1 means that the strategy 
performs as a random guess, and L = means a perfect prediction. If L > 1, 
the strategy is performing exactly the opposite to how it should. Thus, the 
exact reverse strategy should be considered, and this will provide the opposite 
results (/^ = 1 - /„, and f^ = l- f^). 



4.2 Reference strategy 



Of course, the random guessing strategy depicted above is only useful as a 
baseline, but does not serve to provide a real significant forecast. In this subsec- 
tion we describe the simplest meaningful forecasting strategy one can consider 
for any system. This will be called the reference strategy, and any forecasting 
procedure more complex than this should render better results. 

The reference strategy consists simply in declaring an alarm som e time after 

each t arget event, and maintaining it on until the next targ et event ( Newman and Turcottd . 



eacn t arget event, and maintaining it on until tne next targ t 
20021: IVazauez-Prada et al.L booi lOonzalez et"an 1200,^ . 



As a general rule, 

the shorter this time, the bigger fa and the lesser f^. Which time is best, 
then? For the minimalist model, we can look for the number of time steps 
n to use with this strategy for obtaining a smaller L. In a previous paper 
f Vazquez-Prada et al. . 2003[ l we observed that effectively, for each A^, there is 
a n that minimizes L. In Fig. 3, the minimum L that can be obtained with 
this strategy is plotted for A^ between 2 and 20, in the curve labeled "Refer- 
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Unattainable zone 



0,0- 



2 4 6 8 10 12 14 16 18 20 

N 



' Random guessing — o— Clones 1 
Reference —1^ Clones 2 

Ideal — O— Clones 3 



Fig. 3. Loss function (L) obtained with the different forecasting methods used in 
this paper, for various system sizes (2 < < 20). A random guessing strategy 
would render L = 1 for any A^, while L = would mean a perfect prediction. The 
shadowed zone is unattainable for any forecasting strategy used in the minimalist 
model, and the strategy that marks its upper limit is called "Ideal" . The "Reference" 
strategy is based only on the series of the largest earthquakes in the model. The 
three strategies labelled "Clones" are based on the synchronization of models with 
the minimalist model whose largest earthquakes we try to forecast. 

ence". This method does not generate any false alarm, nor take into account 
the occurrence of earthquakes smaller than the characteristic ones. The only 
information required is the statistical distribution fprobability dis tribution 



function) of the duration of the cycles (jVazquez-Prada et al.l . 120031 1 . Taking 



into account the effects of smaller earthquakes, the forecast can be modes tly 
improved in the model J v..n„..-P,.H„ .1 1 h....l., .1 1 TO 



4.3 Ideal strategy 



As the minimalist model is very simple, it is possible to explore its maximum 
predictability. The ideal strategy needed for getting this result, unlike the two 
previously described, is model-specific. It is deduced in Appendix A. This ideal 
result could only be obtained if we could "see" inside the model to check at 
each time step which cells are occupied and which are not. Thus it requires 
a perfect knowledge of the system, and equivalent strategies cannot be used 
with actual faults where we cannot know the detailed state of stress and 
strain. In Appendix A it is deduced that the alarm should be declared at the 
instant in which N —1 cells of the model are full (just at the beginning of the 
plateau with seismic quiescence commented on in Section 2.2). Then, it should 
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be maintained on until the next characteristic earthquake. This is a no-error 
strategy (/e = 0, and L = fa). As the model is stochastic, fa is not zero; a 
minimum alarm time is needed to forecast all the characteristic earthquakes. 
It is given by = L = N/ (n) , where (n) is the average duration of the cycles 
(which depends on A^). This L is also plotted in Fig. 3, in the curve labelled 
"Ideal" . This is the rigorous minimum L that can be obtained in the model. 
A good forecasting strategy should produce a L lower than the "Reference" 
curve and as close as possible to the "Ideal" curve. 



5 Synchronization-based forecasting 



In this section we will describe the novel forecasting method based on the 
synchronization between models, obtained by imposing the rupture area of a 
minimalist model onto other s imilar models. This s ection expands and com- 



plements our previous results (jGonzalez et al.l . 12004^ . 



We will try to forecast the characteristic earthquakes generated by a mini- 
malist model with A^ cells. This model will be called master. We will consider 
this master as if it were an actual fault, from which we can know the rupture 
area of its earthquakes (equivalent to the number of cells broken, k), but not 
the strain or stress at depth (equivalent to the occupation state of the model 
cells). As in an actual fault, we cannot change the state of the master at any 
moment. 

In this forecasting me thod we will use other models, which we call clones 



( Gonzalez et al.L |200J). These are equivalent to the models that a scientist 
devises for forecasting the future evolution of the fault. We will modify their 
evolution at will, and their governing rules will be different than those of the 
master. In this paper, for simplicity, we will consider that the clones are also 
arrays of A^ cells. The average duration of the earthquake cycle in the model 
(average recurre nce interval of the charact eristic earthquakes), (n), strongly 
depends on A^ ( Gomez and Pachecol . 2004 l. Choosing a different A^ for the 



clones will imply a different loading rate of the cells and a different average 
recurrence interval of the characteristic earthquakes in the clones than in the 
master. These effects would require further tuning of the clones, which would 
complicate the following discussion. 

Let us describe in the following paragraphs the general outline of the proce- 
dure. We will use a total of Q clones, that will be loaded (one particle per 
time step and per clone) at the same time as the master, but randomly and 
independently to the master and to each other. We will apply some procedures 
for partially synchronizing the clones with the master. Namely, if in a given 
time step the master does not generate any earthquake, we will oblige the 
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clones not to generate any earthquake either. And if the master does generate 
an earthquake, we will force the clones to reproduce the rupture area of this 
earthquake, as described below in more detail. Note that, although the master 
and the clones are driven simultaneously, the effects of the master are dealt 
with first. 



Why use several clones? The master and the clones are all stochastic, so each 
one evolves with time in a different way. By using several clones, we can take 
into account a broad range of possible evolutions. By using only one clone, 
we could not be very sure that it is satisfactorily mimicking the evolution of 
the master. However, if several of these Q clones are in the same state, then 
it is more probable that the master is also in that state. If the clones were 
deterministic, only one would be required. 

We have commented before (Section 4.3 and Appendix A) that the ideal fore- 
casting strategy for the minimalist model will be to declare the alarm just 
when — 1 cells of the model become occupied. Then the master enters the 
stage of seismic quiescence, or plateau, and the next earthquake is the char- 
acteristic one. We will try to determine this ideal instant as well as possible 
with the clones. For this, we will use a "democratic" procedure: we will declare 
an alarm when a minimum of q clones "vote" (become occupied to a certain 
threshold, described below). Later on we will explore the combinations of Q 
and q that render the best results. Once the alarm is declared, it is maintained 
on until the next earthquake in the master. If it is a characteristic one, this is a 
successful prediction. Its rupture is imposed on the clones (so we reset all the 
cells of the clones to empty) and a new cycle starts. If the next earthquake is 
not a characteristic one, this represents a false alarm. We will disconnect the 
alarm, and impose the rupture on the clones as is done with any other earth- 
quake. Of course, if a characteristic earthquake takes place when the clones 
have still not declared the alarm (when less than q clones have voted), this is 
a prediction failure. If the clones declare an alarm in the same time step in 
which the master generates a characteristic earthquake, we also consider this 
as a prediction failure. 

The exact rules for driving the clones will follow one of the three approaches 
commented on below. Each approach imply a different knowledge of how the 
master works, and a different way of imposing the rupture area on the clones. 
They are depicted in Fig. 4 and described as follows: 



1^1) This first approach will indicate which is the best result that can be 
obtained with the synchronization-based forecasting. For t his reason, the 



clones are indeed minimalist models identical to the master (Gonzalez et al 



20041 ) . The clones are loaded only if the master does not generate an earth- 



quake in that time step. We know that in this case the particle in the 
master has gone to one of the cells i > 2, so the particles in the clones 
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Fig. 4. Sketch that shows how the rupture area of an earthquake in the master 
model (M) is imposed on the clones (C), for each of the three synchronization-based 
approaches. In this example, the master generates an earthquake with rupture area 
A; = 3. In the first approach (CI), the first k + 1 cells rupture and will be reset to 
empty. In the second one (C2), this occurs only with the first k cells. In the third 
approach (C3), this happens to k occupied cells chosen randomly. The first cell of 
the clones can be occupied only in the second and third approaches. 



will be randomly thrown to the cells i > 2. We also consider as known 
that, just after an earthquake with rupture area k, the first k + 1 cells in 
the master, for sure, are stable (the k just broken plus the one that acts 
as a barrier for the rupture). Thus, if the master generates an earthquake 
of size k, we will reset to empty the first k + 1 cells of the clones. A clone 
votes when — 1 of its cells are full. 

(2) In this second approach, we are more ignorant about how the master 
works. At every time step we will throw the stress particles to any of the 
cells in the clones. If the master generates an earthquake of size k, we 
only know which cells have ruptured, so we will reset to empty only the 
first k cells of the clones. A clone votes when its A^ cells are full. 

(3) In the third approach we know even less. At every step we will throw the 
stress particles to any of the cells in the clones. When an earthquake takes 
place in the master, we only know its size, and thus its rupture area, k, 
but not exactly which cells have ruptured. Thus, we will randomly empty 
k occupied cells of each clone. If the clone has less than k occupied cells, 
all are emptied. A clone votes when its A^ cells are full. In this approach 
the positions of the cells in the clone are irrelevant. Each c lone is thus 
equivalent to the so-called box model ( Gonzalez et al. . 2005[ l. 



Note that, ideally, the clones should have the same number of occupied cells as 
the master. For this reason, as a way to measure the degree of synchronization 
between a clone and the master, we used the fraction of time, r, during which 
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Fig. 5. (a) The fraction of time, r, during which two models have the same number 
of occupied cells. This depends on whether they are two independent masters (tq) 
or a master and a clone (governed by one of the three different synchronization 
approaches: ri, T2 and rs). (b) The relative improvement, defined as Tu/tq. 



both of them have the same number of occupied cells ({Gonzalez et al.L |2004 ) 



If two independent masters run simultaneously, they have the same number 
of occupied cells, just by chance, during a certain r. When a clone and a 
master are compared, this r greatly increases, as shown in Fig. 5: partial 
synchronization is achieved. The best results, as expected, are achieved with 
the first of the three approaches. 



The results of /a, /e, // and L = fa + fe, for different values of Q and q can 
be plotted as in the diagrams of Fig. 6. In this figure we have plotted only 
results corresponding to the first of the three approaches and = 20, but 
similar figures, with the same overall properties, can be drawn for the other 
two approaches and for any (see below). There are simple trends in these 
graphs. In Section 3 we noted that, in general, a forecasting strategy that 
produces lower /g tends to produce higher fa and //. If Q is fixed (same row), 
the greater the q, the later the alarm is declared, so and // are lesser and 
/e is greater. If q is fixed (same column), the greater the Q, the earlier the 
alarm is declared, resulting in the opposite trend. 

We are interested in finding the combinations of Q and q that minimize L. 
The interesting fact is that the sum fa + fe shows a rectilinear "valley" for 
certain combinations of Q and q, marked with squares in the graph of Fig. 6. 
This valley goes down as Q and q increase. In Fig. 7 it can be observed that 
the valley goes down indefinitely, tending to a lowest asymptotic value of L. 
We estimate this value, as a function of Q, with a three-parameter exponential 
fit of the form F = a exp[b/{Q + c)], where a, b, and c are parameters. The 
value of a is the asymptotic one for Q oo. This value is represented, for 
each N, in Fig. 3. The fa, ff and fe also have asymptotic trends along this 
valley of L, also plotted in Fig 7. They can also be fitted with the same kind of 
three-parameter distribution, to estimate their asymptotic values as Q oo. 
A nice property is that, as shown in Fig. 7 for a certain case, these forecasting 
approaches predict most of the characteristic earthquakes (/e is low), and have 
a very small fraction of false alarms. Note also that a few tens of clones already 
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Fig. 6. Fraction of alarm time (/a), fraction of false alarm time (//), fraction of 
errors (/e), and loss function (L) obtained with the the first synchronization-based 
forecasting approach, for A'^ = 20 and different numbers of clones (Q) and votes {q). 
The squared cells mark a rectilinear valley in the values of L. 




50 100 150 200 

Q, number of clones used 



250 



Fig. 7. Loss function (L), fraction of alarm time (/a), fraction of false alarm time (fj) 
and fraction of errors (/e) obtained with the first synchronization-based forecasting 
approach, for = 20 and different numbers of clones (Q) along the rectilinear 
valley observed in L in Fig. 6 (the first five points of each curve correspond to the 
cells marked in that figure). 

render results close to the asymptotic ones. 



As can be noted in Fig. 3, the synchronization-based strategies perform much 
better than a random guess, and also much better than the reference strat- 
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Fig. 8. Loss function (L) obtained with the three synchronization-based forecasting 
approaches, for N = 100 and different numbers of clones (Q) and votes {q). The 
squared cehs mark rectihnear vaheys in the values of L. 



egy described in Section 4.2. Their results are intermediate between the ideal 
forecast and the reference one. The second and third synchronization-based 
approaches give only slightly greater L than the first one. The differences 
are large only for small A^. Although the first approach synchronizes more 
efficiently each individual clone with the master (Fig. 5), this effect is com- 
pensated by using many clones. 



To assess the performance of the method with larger systems, we plot in Fig. 8 
the results of L for the three approaches, for = 100 and up to 60 clones. As 
occurred for smaller A^, a rectilinear valley is observed in the graphs, and this 
tendency can be extrapolated to estimate the asymptotic value of L. Note that 
the results for the first and second approaches are almost identical (although 
L is slightly larger in the second approach). With the third approach, L decays 
to its asymptotic value more slowly (the valley floor has a smaller slope, so 
more clones are needed to achieve a given low L). The asymptotic values of 
L, however, are very similar in the three cases (0.298 for the first and second 
approaches; 0.306 for the third one). Note that these values are smaller than 
for A^ = 20, as expected from the trend observed in Fig. 3. The fa, fe and // 
show trends similar to the ones described for Figs. 6 and 7. The asymptotic fe 
is very low (0.062, 0.075 and 0.071 for the first, second, and third approach, 
respectively). 



Another way to measure the synchronization of the clones with the master is 
drawn in Fig. 9 for A^ = 20. The ideal strategy (Section 4.3 and Appendix A), 
would be to declare the alarm just when A^ — 1 cells of the master are full. 
The figure shows how a single clone declares the alarm around that moment, 
but a group of clones does a much better job. 
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Fig. 9. Probability that a clone or thirty clones (with q = 5; the uppermost squared 
cell in L in Fig. 6) declare an alarm in a given time, around the instant when it 
should for obtaining the ideal forecast. 



6 Discussion and conclusions 



In this paper we have tried to provide some insight into how to synchronize 
numerical models with seismic faults, in order to better forecast large earth- 
quakes in them. The idea is that, although we can rarely measure the stress 
and strain in actual faults, we can estimate the rupture area and coseismic 
displacement of their earthquakes. If we force a calibrated model to reproduce 
every earthquake rupture of the fault it simulates, probably the model will 
be synchronized with the fault. Then it could be used to forecast the future 
evolution of the fault, includin g future large earthquakes. This idea is not com- 
pletely new: e.g. IwardI (|2nnnh forced a model to reproduce a large-earthquake 
rupture and run the model forward to see what could happen iri the fu ture. 



The results of this paper expand on earlier ones ([Gonzalez et al.L 12004 ). and 



are still only theoretical, but fully quantitative. We demonstrate that it is 
possible to partially synchronize numerical fault models between themselves, 
and use this to forecast synthetic earthquakes. 



One of the models, called the master, evolves freely. We consider it as an actual 
fault, from which we can know the rupture area of its earthquakes, but not 
the strain or stress at depth. Our goal is to forecast the largest earthquakes it 
generates. In the synchronization-based forecasting, we use several other mod- 
els, called clones, similar to the master (calibrated to have the same average 
recurrence interval of large earthquakes that the master has). These clones are 
equivalent to the models that can be devised to simulate a seismic fault. They 
are run simultaneously and independently to the master and to each other. We 
force them to reproduce the series of earthquake ruptures of the master, and 
this makes them partially synchronized with it. In simple words, if the master 
does not generate an earthquake, we preclude any earthquake in the clones; 
if the master does generate an earthquake, we impose the same rupture area 
on the clones. When several of the clones indicate that a large earthquake is 
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impending in the master, we declare an alarm. This efficiently predicts most 
of the largest earthquakes of the master, with a relatively low fraction of total 
alarm time and few false alarms. These results are robust: they are almost 
the same when the exact rules for imposing the earthquake ruptures vary, and 
this good performance is observed along the whole range of model sizes con- 
sidered. This synchronization-based forecasting o utperforms other procedures 



based only on the ear thquake series of the model ( Vazquez- Prada et al. . 20031 : 



Gonzalez et al.L l2004^ 



The master and the clones are stochastic (random), so each individual clone 
is only partially synchronized with the master. However, when several clones 
are in the same state, then it is more likely that the master is also in this 
state, so the group of clones makes a much better forecasting job than only 
one clone does. If the clones were deterministic, as a general rule only one 
would be needed; more clones would have identical evolutions if run with the 
same initial conditions. 



The procedure developed here is a kind of ensemble forecasting, in which 
several models are run to obtain a better picture of how a system wi ll evolve. 



This concept is used in atmospheric forecasting ( Palmer et al. . 2005[ ): several 



models are run simultaneously, and their average result has a larger forecastin g 
ability than that of an individual model (Houghton, 20021 : Palmer et al. . 2005| ). 
Each model in this approach has slightly different initial conditions, to take 
into account measurement errors and then to represent one possible state of the 
atmosphere, among various possibilities. In our approach, each clone marks 
a possible state of the master among a range of possible options. Several 
deterministic clones could also be used with different initial conditions. 



Our procedure also shares some si r nilari t ies with certain earthquake forec ast- 
ing algorithms ( Kossobokov et al. . 19991 : Keilis-Borok and Soloviev . 2003[ ). in 
which several seismicity functions are evaluated in real time. When several 
of these functions indicate that a large earthquake is probable, an alarm is 
declared. In our approach, the clones are performing a role similar to these 
functions, monitoring what is happening in the master. 

Our proposal is that a possible way to synchronize more complex, calibrated 
models with real faults might be to force them to reproduce the past series 
of earthquakes (with the same rupture area and/or coseismic displacement). 
This would need to be tested in the future. Also it will be possible to test 
whether this procedure works in the forecasting of synthetic earthquakes in 
other models. 



Forci ng the models to reproduce only one large observed rupture (as in [Ward , 
2000l ) probably is not enough (this is certainly the case in our stochastic 
model). Complex and chaotic systems, such as the lithosphere, are very sensi- 



20 



tive to initial conditions. Forcing the model to reproduce only one rupture is a 
necessary and laudable first step, but probably does not constrain the initial 
conditions sufficiently. We propose that every observed rupture, albeit small, 
should be considered. Small earthquakes are much more frequent than large 
ones, thus providing much more data. Moreover, they provide insight into 
the mechanical state of the crust ( Seeber and Armbrusteil . 2000( ) and into the 
mechanics of earthquake rupture ()Rubin . 2002| ). Their location may indicate 
the patch of the fault plane which i s experiencing higher stres s es an d is likely 
to rupture in the next large shock ( Schorlemmer and Wiemer . 2005[ ). Finally, 
they are important in the transfer of stress w ithin the lithosphere, and in 
earthquake triggering (|Helmstetter et aLllimih . 
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A Deduction of the ideal forecasting strategy 



In this appendix we will deduce the ideal strategy outlined in Section 4.3. 
This strategy renders the lowest (best) value oi L = fa + fe achievable in the 
minimalist model. 

For this reasoning we would consider every cycle of the model as composed 
of two independent and consecutive stages. The ffist, that will be called the 
loading stage, starts just after the occurrence of a characteristic earthquake. 
During this stage the total number of occupied cells grows, but not in a mono- 
tonic way, because the particles may land in already occupied cells (and then 
be dissipated), and also because of the occurrence of non-characteristic earth- 
quakes (Fig. 2). When — 1 cells (all but the ffist one) become occupied, 
this ffist stage ends and the second stage, that will be called the hitting stage 
(or plateau in the occupation), starts. In this second stage, the system resides 
statically in the state of maximum occupancy (Fig. 2) until a particle arrives 
at the ffist cell. Then, a characteristic event occurs, all the cells are emptied, 
and a new cycle begins. The hitting stage can be mathematically treated as a 
form of Russian roulette. 
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Both the time spent by the system in the loading stage, x, and in the hitting 
stage, y, are statistically distributed. The distribution of y, denoted by Ad/), 
is geometric. Considering that, in each time step, the probability of hitting 
the first cell is p = and its complementary is g = 1 — it follows that 

1 / 1 \^"^ 

^2(y) = T7 1-T7 , (A.l) 



whose mean is 

{y) = N, (A.2) 
and whose standard deviation is 

(A.3) 



The time elapsed between consecutive characteristic events has been denoted 

by n, which is st atistically distrib uted according to the function PAr(n) (|Vazquez-Prada et al, 




200^ 



J, 2003; .Gomez and Pachecol . Eo04 ) . Because the variables x and y are in- 



dependent, the mean length of the cycles (n) is the sum of the mean lengths 
of the two stages: 

(n) = {x) + {y). (A.4) 



It is clear that the best L would be obtained only if we knew the state of 
occupation of the system and could mark, for each cycle, the instant at which 
the stage of loading concludes. In this case, fe = 0, but because the hitting 
stage is completely stochastic, fa (and thus L) cannot be nil. 

Let us explore the result of L obtained if we turn the alarm on at a given 
value y = yo within the second stage of all the cycles. With this strategy, the 
fraction of errors is given by 

yo 

/eM=E^2(2/), (A.5) 



and inserting Eq. A.l we obtain 

/e(2/o) = l-g^°. (A.6) 
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With respect to to the fraction of alarm, its form is 

oo 

faiyo) = ^° , . , (A.7) 

(x) + {y) 

and inserting Eq. A.l, we get 

Um)-j^- (A.8) 

Note the important contribution of the first stage of the process in the de- 
nominator. Thus, the specific form of the loss function is 



It is noteworthy that in the absence of the first stage, i.e. in the hypothesis 
of a pure geometric distribution, the value of L would be 1, not dependent on 
the value of t/q. In this sense, the geometrical and the Poisson distributions are 
equivalent. The minimum value of L in Eq. A. 9 as a function of yo is obtained 
for yo = 0, i.e. just after the end of the first stage, when the — 1 upper cells 
of the system are full. And this minimum value is 

N 

Lmin = — r- (A. 10) 

(n) 



This result constitutes a rigorous lower bound for the expected accuracy of 
any forecasting strat egy in the minimalist niodel . For this model, (n) increases 
rapidly as N grows ( Gomez and Pachecol . 2004 ). This implies that the min- 
imum L, obtained with this optimal forecasting strategy, decreases as N in- 
creases, as shown by the curve labeled as "Ideal" in Fig. 3. That is to say, 
minimalist models with more cells are more predictable. This is consistent 
with the fact that the time series of characte ristic earthquakes is more peri- 
odic for larger ( Gomez and Pachecol . 2004). 
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